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We extend the generalized Langevin model, ^ originally developed for the Lagrangian fluid particle 
velocity in constant-density shear-driven turbulence, to variable-density (VD) pressure-gradient- 
driven flows. VD effects due to non-uniform mass concentrations (e.g. mixing of different species) 
are considered. In the extended model large density fluctuations leading to large differential fluid 
accelerations are accounted for. This is an essential ingredient to represent the strong coupling 
between the density and velocity flelds in VD hydrodynamics driven by active scalar mixing. The 
small scale anisotropy, a fundamentally "non-Kolmogorovian" feature of pressure-gradient-driven 
flows, is captured by a tensorial stochastic diffusion term. The extension is so constructed that it 
reduces to the original Langevin model in the limit of constant density. 

We show that coupling a Lagrangian mass-density particle model to the proposed extended ve- 
locity equation results in a statistical representation of VD turbulence that has important benefits. 
Namely, the effects of the mass flux and the specific volume, both essential in the prediction of VD 
flows, are retained in closed form and require no explicit closure assumptions. 

The paper seeks to describe a theoretical framework necessary for subsequent applications. We 
derive the rigorous mathematical consequences of assuming a particular functional form of the 
stochastic momentum equation coupled to the stochastic density field in variable-density flows. Our 
aim is to develop a joint model for variable-density pressure-gradient-driven turbulence and mixing, 
such as occurs due to the Rayleigh- Taylor instability. A previous article^ discussed VD mixing and 
developed a stochastic Lagrangian model equation for the mass-density. Second in the series, this 
article develops the momentum equation for VD hydrodynamics. A third, forthcoming paper will 
combine these ideas on mixing and hydrodynamics into a comprehensive framework: it will specify 
a joint model for the coupled problem and validate it by numerically computing joint statistics of a 
Rayleigh- Taylor flow at several Atwood numbers. 

Keywords: Probability density function method; Langevin equation; Variable-density turbulence; Pressure- 
gradient-driven flows; Small-scale anisotropy 
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I. INTRODUCTION 

Haworth and Popc^ developed the stochastic model for 
the Lagrangian fluid particle velocity, w*, in constant- 
density shear-driven turbulence, 



du* 



/godt + G.y {v* - V,) dt + {G^ef^ AW,, (1) 



where Qo, p and e are the constant density, pressure and 
the dissipation rate of turbulent kinetic energy, respec- 
tively. The coefficients, Gij and Co, are specified by the 
particular model of choice. Eq. (1) is a stochastic differ- 
ential equation (SDE) of the diffusion type, where dWi{t) 
is an isotropic Wiener process'^ and the overbar denotes 
the ensemble average. 

Eq. (1) is used to represent the probability density 
function (PDF) of velocity. In PDF methods'* a trans- 
port equation is solved for the joint PDF of several flow 
variables. In most cases the equation is of Fokker-Planck- 
type and its numerical solution amounts to following a 
large number of Lagrangian particles in a Monte-Carlo 
fashion, whose properties are governed by an equivalent 
system of SDEs, such as Eq. (1) for the particle veloc- 
ity. These computational particles do not correspond to 
real fluid particles, but the statistics of an ensemble are 
representative of the local statistical behavior of the flow 



over many realizations. 



A. Variable-density pressure-gradient-driven 
turbulence 



Our goal is to extend Eq. (1) to variable- density (VD) 
pressure-gradient-driven turbulence (PGDT), a class of 
flows that is fundamentally different from constant- 
density shear-driven flows. A classical example of PGDT 
is the Rayleigh- Taylor (RT) instability of an interface 
between two fluids of different densities in a gravita- 
tional field that points to the opposite of the density 
gradient. RT flows are important in several geophysi- 
cal, astrophysical and engineering applications. 

Boussinesq limit. If the densities of the mixing flu- 
ids are commensurate, PGDT is treated in the Boussi- 
nesq limit: the density fluctuations are small compared 
to the mean density and are only important in the body 
force term of the Navier-Stokes equation. A discussion of 
RT-turbulence in the Boussinesq case, with a self-similar 
analysis at the second moment level, is given by Ristor- 
celli and Clark. ^ 

Variable-density case. In situations where the fluids 
have vastly different densities the term variable-density is 



2 



LA-UR 10-05227, vO.9, Accepted in Journal of Turbulence, Jan 7, 2011 



used to distinguish it from the Boussinesq case. A com- 
monly used measure of VD effects is the Atwood number, 

A^^^ ^ ^ = i±4^ (2) 
Q2 + Qi Qi I- A 

where qi < Q2 denote the constant densities of the pure 
fluids in a binary mixture. ^ ss corresponds to the 
Boussinesq case, while ^4 — >■ 1 to largely disparate species 
densities, Qi <C Q2- Mixing at high Atwood numbers, 
compared to the Boussinesq case, is accompanied by sev- 
eral new effects: (1) the advection term in the Navier- 
Stokes equation gives rise to cubic non-linearities, (2) ad- 
ditional non-linearities in the molecular diffusion terms, 
(3) a dynamic mean pressure gradient, and (4) the asym- 
metry of the mixing layer. In these flows the mean spe- 
cific volume becomes an important independent variable. 
As a consequence of these new phenomena both hydro- 
dynamics and mixing exhibit qualitatively different be- 
havior compared to the Boussinesq case.^"^"^ 

Challenges associated with VD-PGDT. The 
above description highlights the tight coupling of the 
density and velocity fields in VD fiows. This poses 
unique challenges, compared to classical (and more ex- 
tensively researched) constant-density shear-driven tur- 
bulence. They can be enumerated as follows: 

1. Mixing- driven: The hydrodynamics of VD-PGDT 
may primarily be driven by material mixing. Be- 
side mean deformation a dominant production 
mechanism of turbulent kinetic energy is due to the 
joint effects of non- uniform fluid concentrations and 
pressure gradients. The (stochastic) density field is 
an active scalar in the momentum equation. 

2. Turbulence asymmetry: At high Atwood numbers 
the mixing process becomes asymmetric. Fluid 
particles representing different instantaneous den- 
sity (and thus different inertia) respond very dif- 
ferently to pressure gradients. As a consequence, 
in a homogeneous Ray leigh- Taylor flow an initially 
symmetric density distribution quickly develops a 
sizeable skewness.^'^^ Similarly, the non-Boussinesq 
behavior of the pressure leads to asymmetry in 
the proflles of the mean velocity and the turbulent 
kinetic energy across inhomogeneous RT mixing 
layers. ^■^ In contrast the Boussinesq case remains 
symmetric. 

3. Non- equilibrium: The flow evolution is highly non- 
equilibrium as the production-to-dissipation ratio 
ranges from hundreds to nearly zero.^ 

4. Transitional: The Rayleigh- Taylor mixing layer 
transitions to fully developed turbulence starting 
from a quiescent state. 

5. Non- stationary: At no point in time can the flow 
be deemed statistically stationary. 



6. Anisotropic: Due to the large external acceleration 
force, anisotropy is important at both large and 
small scales. ^"'^ 

The fundamental physics of the flow evolution must be 
captured at several different stages. As an example, in an 
inhomogeneous Rayleigh- Taylor flow,^^'^^ two pure fluids 
are initially separated by a thin perturbed interface in a 
quiescent state: the heavy fluid lies on top of the light 
one with gravity acting downwards. Due to the unsta- 
ble configuration, bubbles and spikes grow and penetrate 
into each other in a laminar stage, which then break up 
into smaller Kelvin-Helmholtz-like eddies. The fiow tran- 
sitions into a fully developed turbulent mixing layer with 
its width growing as long as pure fiuid is entrained at the 
edges. Another example is the homogeneous Rayleigh- 
Taylor mixing layer, ^'^"^ in which the process starts from 
a quiescent state with random blobs of two pure fluids (in 
a closed computational box), corresponding to the center 
of an inhomogeneous RT layer. Potential energy is grad- 
ually converted to kinetic energy during transition to a 
fully developed turbulent state. After reaching a peak in 
kinetic energy, a viscous decay due to small scale dissi- 
pation ensues. 

Modeling challenges of VD turbulence. From 
the viewpoint of statistical turbulence modeling the main 
challenges are: 

• To devise a mutually consistent representation for 
hydrodynamics and mixing. The model must cor- 
rectly take into account the effect of the active 
scalar, the density, on the velocity field, and the 
mixing of the scalar by the turbulent velocity field. 

• To represent very different flow- evolutionary stages 
within the same method. The model must evolve 
the flow from a quiescent initial condition through 
an initially laminar stage, via the highly non-linear 
process of transition to fully developed turbulence, 
followed by a possible decay if pure fiuids are no 
longer entrained. 

• To provide a higher level statistical description than 
current VD moment closures. We require the 
model to predict the full one-point one-time prob- 
ability distribution function of the active mixing 
field (the density), due to its fundamental role in 
mixing-driven VD flows. 

• To formulate a simple method that captures the es- 
sentials. The numerical method must be computa- 
tionally inexpensive, so that it can be used as an 
engineering closure. 

PDF methods seek to compute the one-point joint 
PDF of the fluctuating flow variables. Reviews on PDF 
methods for turbulent flows are provided by Pope,"* 
KoUmann,^^ Dopazo^^ and by Haworth.^'' The most 
widely used velocity PDF models are variants of Eq. (1): 
by specifying Gij and Co differently several models have 
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been constructed (see e.g. Ref. 18) that are routinely used 
today in combustion simulations.^'' To highlight the chal- 
lenge of PDF methods for VD turbulence, the fundamen- 
tal assumptions of the standard PDF methodology for 
hydrodynamics, Eq. (1), are enumerated as follows: 

1. Uniform- density formulation. Eq. (1) has orig- 
inally been developed for constant-density flows, 
go = const. In low-speed combustion the fluid den- 
sity is treated as a function of the species concentra- 
tions and temperature, g — g{y,T), and the mean 
density, p, is used in place of go. This non-constant, 
but from the viewpoint of the momentum equation 
still uniform-density case, does not take direct ef- 
fects of a stochastic variable density field on the 
particle momentum into account. This is most eas- 
ily seen in the pressure-gradient term of Eq. (1) and 
less directly in the other two terms which will be 
discussed in more detail. 

2. Small scale isotropy. With the isotropic stochas- 
tic diffusion term, SijdWj, Eq. (1) adheres to 
Kolmogorov's hypothesis of local isotropy by con- 
struction. Situations where the hypothesis is 
known to fail are non-equilibrium, highly-distorted, 
moderate- i?e or VD flows. The fluctuating velocity 
field in these cases may exhibit significant level of 
anisotropy at the small scales. 

3. Fully developed turbulence. Eq. (1) assumes a fully 
developed, i.e. high-i?e, flow with isotropic incrtial 
range scaling.^ 

The above general discussion makes it clear that devel- 
oping a PDF model for VD flows is a largely unexplored 
territory. Due to the wide variety of challenges associ- 
ated with non-stationary variable-density flows, the de- 
velopment is documented in a series of three papers, as 
described in the following subsection. This article is the 
second in the series. 



B. Objectives of the article 

Our aim is to develop a model that represents the 
temporally evolving joint PDF of density and velocity 
in variable-density pressure-gradient-driven turbulence. 
The development is carried out in three parts: 

1. Development of a VD material mixing model. A 
previous article,^ devoted to the density equation, 
addressed the challenges associated with active 
scalar mixing in VD flows and developed a La- 
grangian stochastic equation to represent conser- 
vation of mass. 

2. Development of the Langevin model for the velocity 
field in VD flows. The extension of the momentum 
equation from constant-density shear-driven flows 
to VD-PGDT is the subject of this paper. Two 



critical ingredients and their rigorous mathematical 
consequences will be discussed: 

(a) The representation of the instantaneous den- 
sity field in the momentum equation, resulting 
in an appearance of the mass flux (and thus 
turbulence production) that does not require 
explicit modeling; and, 

(b) The introduction of a new tensorial stochas- 
tic diffusion term representing small scale 
anisotropy. 

3. Joint PDF model specification and validation. The 
strongly coupled nature of active material mixing 
and hydrodynamics in VD-PGDT will require a 
joint model. This is the subject of a forthcom- 
ing article, which will combine the ideas about 
mixing and hydrodynamics into a comprehensive 
modeling framework. 

C. Navier-Stoltes equation for variable-density flow 

It seems appropriate to specify, at the outset, the dy- 
namic level of hydrodynamical approximation that is to 
be modeled. 

We consider variable-density flows, where differences 
in fluid density arise solely from non-uniform mass con- 
centrations, e.g. due to mixing of different fluids. Den- 
sity variations originating from pressure and temperature 
changes are neglected at this time, g(y,p,T) ~ g{y)- 

For this class of variable-density flows, the Navier- 
Stokes equation, governing the instantaneous fluid par- 
ticle velocity, Vi, in the presence of viscous and external 
forces, is written symbolically as 

dwi = {gi - vp,i -|-/iuV^Wi) dt, (3) 

with the specific volume, 

V = l/g. (4) 

Here gt, p, g and ^ denote the acceleration force per unit 
mass, the pressure, the density and the constant dynamic 
viscosity, respectively. Although in VD flows the viscos- 
ity is non-uniform, at this time we take /i w const, as 
we are interested in exploring the new simpler physics 
(VD turbulence) without the complicating issues of non- 
uniform viscosity. Accordingly, in Eq. (3) the viscous 
force is approximated as v[^{vi^j -l-Uj^i) — 2/3/iWfc,fc(5ij],j « 

The constant-density counterpart of Eq. (3), to which 
the Langevin equation (1) is a model in turbulent flows, 
reads 

dui = [gi - p,i /go + ^/ goV'^Vi) dt. (5) 

In the following we will denote a constant density by go 
and the variable one by g. 
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D. Outline of the paper 

To set the stage, the generaUzed Langevin model, 
Eq. (1), developed for shear-driven constant-density 
flows, is reviewed in Sec. II. This is to point out its most 
important underlying assumptions and so that the later 
development can highlight the major differences in VD 
flows. The Langevin equation is extended to VD-PGDT 
in Sec. III. Sec. IV investigates the new features of the 
variable-density SDE by deriving its ensemble moment 
equations and comparing them to those derived from the 
Navier-Stokes equation. Sec. V does the same for the 
Favre moments. This is to put the proposed PDF for- 
mulation into context with Favre VD moment closures. 
Sec. VI discusses the VD PDF model in the Boussincsq 
limit. The important characteristics of the formulation 
and the main results are summarized in Sec. VII. 



II. REVIEW OF THE LANGEVIN MODEL FOR 
CONSTANT-DENSITY SHEAR FLOWS 

This section reviews the generalized Langevin model 
and its underlying assumptions. 

Decomposing the velocity and pressure into mean 
and fluctuating parts, the Navier-Stokes equation for 
constant-density flows, Eq. (5), takes the form 



dwj = (g,, - p,j /Qf) + ^/ qqV Vi) dt 



(6) 



with the fluctuation defined as y' ^ y — y. 

A widely used model equation for the particle veloc- 
ity, V* , in constant-density shear flows is the generalized 
Langevin model (GLM),^ 



dv* = {gi - p,i / go + ^/fi'oV^Wi) dt 

+ G,j {v* - V,) dt + {Coef^ dW,. 



Strictly speaking, turbulence is not a Markovian 
process, but the fluctuations in the inertial sub- 
range, for which Eq. (7) has been developed, can 
be closely approximated by Markov processes. 

2. Linear dependence: through the linear drift term 
Gij{v* — Vj)dt the effect of fluctuations of the sur- 
rounding fluid is modeled as a linear function of 
the velocity. This is inconsistent with the quadratic 
dependence of the fluctuating pressure gradient on 
the fluctuating velocity. However, linearity (as an 
approximation) may be justified by (1) the correct 
behavior of the equation in homogeneous shear tur- 
bulence, i.e. an arbitrary velocity PDF relaxes to 
a joint normal, and (2) that realizability of the 
Reynolds stress tensor is automatically satisfied 
provided that Co is non-negative and Gij and Cq 
arc bounded. ^'^'^ 

3. Local isotropy: the isotropy of the diffusion term 
{Coey^^SijdWj implies isotropy at the small scales, 
thereby ensuring consistency with Kolmogorov's 
hypothesis of local isotropy. 

Since the first drift term of Eq. (7) represents the effects 
of the mean forces on the particle and the last two terms 
do not affect the mean, the equation is consistent with 
the mean of the Navier-Stokes equation (5).^^'^^ A direct 
consequence of the linearity of the second drift term in 
V* and the independence of v* of the stochastic term is 
Gaussianity of the joint velocity PDF.'^ 

Comparing Eqs. (6) and (7) shows that the last two 
terms in Eq. (7) jointly model the combined effect of the 
fiuctuating pressure gradient and viscous dissipation: 



-p'ri /qo + A^/e'oV dt 



G,,{v* ^Vj)dt + iCoef^dW,. 



(9) 



(7) 



An important consistency condition on the coefficients 
Gij and Co is^ 



In the following, the star superscript (*) will denote a 
model for an instantaneous quantity, such as v* for Vi. 

In Eq. (7) Gij is a second-order tensor which, in 
shear fiows, is assumed to depend on local values of the 
Reynolds stress, v[v'j, the dissipation rate of turbulent ki- 
netic energy, e, and the mean velocity gradient, Vij . Co is 
a positive constant and dWi (t) is a vector- valued Wiener 
process^ with zero mean and covariance dH/idWj = dtSij. 
The statistics on the right hand side of Eq. (7) are under- 
stood to be evaluated at the particle position, x*, gov- 
erned by 

dx* ^ v*dt. (8) 
Three fundamental assumptions undcrly Eq. (7): 

1. The Markov property,^ which assumes that the ve- 
locity can be described by a Fokkcr-Planck equa- 
tion and, equivalently, by a SDE, such as Eq. (7). 



-Co 



0, 



(10) 



which ensures that no spurious turbulent kinetic energy, 
k = v[v[/2, is created in homogeneous turbulence by the 
model. Eq. (10) is obtained from comparing the evolution 
equations that govern k in homogeneous turbulence: the 
one derived from the Langevin model, Eq. (7), and the 
one from the Navier-Stokes equation (5). The constraint 
in Eq. (10) ensures the same form for the k equation: 



dk 
'dt 



GLM 



dk 
'dt 



(11) 



NS 



where V is the shear production. The simplest way to 
define Gij to satisfy Eq. (10) is then 



Gy - - ( - + -Co ) -S^j, 



(12) 
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resulting in the simplified Langevin model (SLM),-"^ which 
corresponds (at the Reynolds stress level) to Rotta's 
model of return-to-isotropy.^^ 

The most important characteristics of the generalized 
Langevin model have been summarized: (1) the model 
equation is consistent with the mean Navier-Stokes equa- 
tion; (2) the mathematical form accommodates Kol- 
mogorov's hypothesis of local isotropy; (3) a realizable 
Reynolds stress model is ensured, provided the coefH- 
cients satisfy certain mild conditions; and (4) the pre- 
dicted joint velocity distribution is Gaussian. 

III. A LANGEVIN MODEL FOR 
VARIABLE-DENSITY 
PRESSURE-GRADIENT-DRIVEN FLOWS 

This section extends the Langevin model, Eq. (7), to 
VD-PGDT and discusses its main ingredients. 

We propose to model the Lagrangian velocity incre- 
ment in YD turbulence by 

dv* = {g, - p„ /g* + iu/q*W'^v,) dt + {v* - {v^)) dt 
+ (0.£)'/' dW, + iM'^^ h,jdW;, (13) 

with positive bounded tensorial diffusion hij, bounded 
functions (j), > and > 0, and independent Wiener 
processes dWi and dWj. 

The SDE (13) differs from (7) in three distinct ways, 
subsequently discussed in more detail in the following 
three subsections: 

A. Instantaneous density. Eq. (13) is intended to be 
coupled to a stochastic equation governing the in- 
stantaneous density field, g*, representing conser- 
vation of mass. 

B. Relaxation to Favre- averaged velocity. The lin- 
ear relaxation term, Gij{v* — {vj))dt, involves the 
Favre-averaged velocity, (vi) = gvi /'g. 

C. Small seale anisotropy. A new tensorial diffusion 
term, {(j)o£)^^'^hijdWj, is introduced. The con- 
stant, Co, in the original isotropic diffusion term 
is exchanged to which together with (jj^ will be 
specified later. 

A. Instantaneous density 

In YD turbulence large density variations play a ma- 
jor role and, compared to the Boussinesq case, the den- 
sity fluctuations can no longer be neglected in the iner- 
tia terms of the Navier-Stokes equation. Extending the 
SDE (7) to YD flows requires a representation of the 
variable density field. Eq. (13) accomplishes this in the 
most profitable way allowed by one-point PDF methods: 
a stochastic density equation, such as discussed in Rcf. 
2, is coupled at the instantaneous (particle) level. 



Comparing the constant-density and variable-density 
models, Eqs. (7) and (13), we see that in the YD case the 
particle density, g*, divides the mean pressure gradient 
and the large-scale viscous terms, 

-Pn/ga ^ -VnlQ*-, (14) 
^i/goV\H i^i/g*V^Vi, (15) 

which has the following consequences: 

• The full density PDF is coupled to the mean forces 
in the momentum equation. Since the governing 
equation for g* provides the full density PDF, the 
mixing state is represented by including the effects 
of all density moments. As Ref. 21 will demon- 
strate, this is crucial in capturing the asymmetric 
PDF of the fluid density at high Atwood numbers. 
Eqs. (14) and (15) ensure the coupling of the den- 
sity PDF to (i.e. the effects of all its moments on) 
the mean forces. 

• The effects of the mass flux on the Reynolds stress 
appear closed. As will be shown in Sec. lY, Eq. (14) 
is the key to represent the effect of the mass 
flux on the Reynolds stress in closed mathematical 
form. This is crucial in variable-density pressure- 
gradient- driven flows, as the mass flux relates to 
a primary mechanism of kinetic energy production, 
which in moment closures requires additional model 
equations. 

The Navier-Stokes equation for YD flows, Eq. (3), with 
the velocity and pressure decomposed, is 

dui = [gi - vp,i +fivV^Vt) dt + {-vp\i +/iwV^w-) dt, 

(16) 

where the instantaneous specific volume, ij = w + w', is 
not decomposed. Comparing Eqs. (13) and (16) shows 
that the combined effects of the fluctuating pressure gra- 
dient and viscous dissipation (which now arc multiplied 
by the instantaneous specific volume, v = \/ g) are jointly 
modeled as 

(- V,. +l^vV\[) dt = G,, (v* - {v,)) dt 

+ (0.e)'^' dW, + {c^nef'^ h.jdW^. 

In other words, in Eq. (13) the last three terms are a 
model, while the representation of the mean forces (mod- 
ulated by the instantaneous specific volume, v* = l/^?*), 
gi —p,i I g* + ^/ g*'S/'^Vi, is exact. This is a consequence of 
the availability of the one-point density PDF, represented 
here via the instantaneous g* . Although the governing 
equation of g* may contain modeling,^ the coupling of 
g* to the mean forces terms of the momentum SDE (13) 
is mathematically exact. As long as the joint density- 
velocity PDF is valid (whose definition is discussed later) 
and the marginal density PDF ensures conservation of 
mass, the statistics involving the density, such as g'v[, 
are finite, consistent and physically realizable. This is 
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discussed further in Sec. IV. In this paper, we do not 
assume a particular functional form for g* , only its avail- 
ability and validity. Rcf. 2 discusses one choice of the 
density equation, which can be used in conjunction with 
Eq. (13). However, we emphasize that the current de- 
velopment is independent of the functional form of the 
density model. 

B. Relaxation to Favre-averaged velocity 

Comparing the constant-density and variable-density 
SDEs, (7) and (13), shows that the linear relaxation term 
in the VD case involves the Favrc-average of the velocity 
field, {vi) ^ Wi/Q as 

^ G,, {v* - {v,)) . (18) 

The same approach has been taken by Delarue & Pope^^ 
to develop a PDF model for high-speed compressible 
shear flows. In the application of Eq. (7) to turbulent 
combustion relaxation to the Favre average is a standard 
procedure. Sec. IV will show that Eq. (18) is instrumen- 
tal in representing the effects of the fluctuating specific 
volume on the mean velocity, important in mixing flows 
with large density differences. 

C. Small scale anisotropy 

In higli-Rcynolds-number shear flows it is generally as- 
sumed that the small scales become isotropic and in- 
dependent of the (anisotropic) large scales where most 
of the turbulence production takes place. Exceptions 
arc strongly distorted flows^*^ or, as recently shown by 
Livcscu & Ristorcelli,^° Livescu et al.^^'^'^ and Chung 
& PuUin,^^ the Ray leigh- Taylor mixing layer. In RT 
flows the buoyancy force has a significant effect on the 
smallest scales, resulting in small scale anisotropy. This 
is apparent in the Reynolds stress anisotropy, bij{k) ~ 
'^'i'^'jl'^'k^'k ~ ^ijl^ at high wavenumbers k,^'^'^^ and in the 
dissipation rate anisotropy, dij =eij /ekk — Sij /3, Ref. 10. 
The above studies show that this is prevalent in the At- 
wood number range A = 0.04 ~ 0.75. We are interested 
in the full Atwood number range of < ^ < 1. Conse- 
quently, assuming small scale isotropy, dij ~ 0, in such 
pressure-gradient-driven flows is not justified. Account- 
ing for the dissipation rate anisotropy is crucial in pre- 
dicting the correct Reynolds stress tensor, whose bud- 
get is directly affected by dij via the small scale dissipa- 
tion term Sij = 2e{dij + Sij/3), see also Eq. (36). The 
anisotropic behavior of e^, governed by dij, is responsi- 
ble for different dissipation rates of the individual com- 
ponents of the Reynolds stress. 

In the constant-density model SDE (7), the small scales 
are isotropic by construction, which is built into the 
stochastic term, (Coey^^dijdWj, in accord with Kol- 
mogorov's hypothesis. As will be shown in Sec. IV, in the 



absence of shear production, a source of (single-point) 
anisotropy in bij can be represented by an anisotropic 
specification for dj, while both &y(/i)ssO and dij{k)~0 
at the small scales. As PGDT exhibits anisotropy at 
both large and small scales (in both bij and dij), the 
variable-density velocity model, Eq. (13), relaxes the 
third assumption implied by Eq. (7) and abandons the 
consistency with Kolmogorov's hypothesis of small scale 
isotropy. It replaces the isotropic diffusion term by the 
sum of an isotropic and an anisotropic (tensorial) diffu- 
sion term as 

(Cos)'/' dW, (0,e)'/' dW, + {cj)oef^ K^dW'^. 

(19) 

The first term in Eq. (19), proportional to is designed 
to account for the effects of the isotropic part of the 
kinetic energy production/dissipation, while the second 
one, proportional to 4>u and hij, for the effects of its de- 
viatoric part. As will be shown in Sec. IV, Eqs. (28-29), 
the stochastic diffusion terms in Eq. (19) arc source/sink 
terms, depending on the sign of the given Reynolds stress 
component. 

The tensorial term in Eq. (19) is assumed to be pro- 
portional to e = efefc/2. This is a conventional procedure 
in low-Reynolds-number flows, such as near walls, ^''■^^ 
and justified here by dijit) being non-zero through all 
the Ray leigh- Taylor flow evolution, see Fig. 16 in Ref. 
10. In other words, the scalar kinetic energy dissipation 
rate, e, is made anisotropic by the tensor hij in the VD 
model. This ensures that the small scales are anisotropic: 
the Lagrangian velocity structure function of the process 
governed by Eq. (13) becomes'^^ 

A, {s) EE [vnt+^'^)-v*mv*it+s)-v*{t)] ^^^^ 

= ((/)/% + 4>ohtkhkj)es, 

with s <^ T <^T , where r and T = \\vk,k\\ ^ denote the 
time scales of the dissipation and the mean deformation, 
respectively. As will be shown in Sec. IV, hij, as a contri- 
bution of the small scales, provides an additional source 
(beside dj) in the budget of the one-point Reynolds 
stress anisotropy, bij. 

Constraints. Up to this point, the functional forms 
of the diffusion terms for VD flows, given by Eq. (19), are 
formulated based on physical insight and mathematical 
consistency. These are: 

• Small scale anisotropy. The anisotropy of RT flows, 
indicated by both bij ^ and dij 7^ at the small 
scales. 

1 /2 

• Correct unit. As 0/ and ipD hij are assumed to be 
non-dimensional, each component of the diffusion 
terms has the unit as that of V edt. 

• Second-order tensor coefficient. The simplest way 
to introduce directional dependence into the joint 
statistics of the fluctuating velocity components is 
to do so via the product of a second-order tensor 
and a vector-valued Wiener process. 
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The cocfRcients 



br, and hij must be specified based 



on the following considerations: 

1. Consistency with the SDE. The coefficient tensor 
should be a square root of a symmetric and non- 
negative semi-definite tensor, Hij = hikhjk- This 
is required mathematically for a SDE, such as 
Eq. (13); to represent a diffusion. ■^^ 

2. Reflect the physical source of small scale anisotropy. 
The diffusion coefficient tensor, hij, should be 
a function of the given source of small scale 
anisotropy. In buoyantly driven flows this may be 
the buoyancy force, g^, in other types of externally 
accelerated or strongly distorted flows, the body 
force, the pressure-gradient force or the mean strain 
rates, responsible for the anisotropy. This ensures 
that the model directly represents the effect of the 
source of anisotropy on the small scales. 

3. Correct distribution. None of hij, (j), or (j)^ should 
be an explicit function of the particle velocity, v* , 
if the Gaussianity of the joint velocity PDF is to 
be preserved. In constant-density homogeneous 
shear flows the velocity distribution is very close to 
a joint normal. "^^ Low-Atwood-number direct nu- 
merical simulations (DNS) show that joint Gaus- 
sianity of the velocity components is also a good 
approximation in RT flows. '^'^ As the drift of the 
Langevin equation (7) is linear and its diffusion is 
not an explicit function of v* , it is an Ornstein- 
Uhlenbeck process,"^ whose statistically stationary 
solution is a Gaussian. As an asymptotic require- 
ment of the VD model Eq. (13), that is to reduce 
to the constant-density Langevin model, the new 
diffusion terms should approach this attribute as a 
limit when g' 0. 

4. Asymptotic small scale isotropy. As an extension 
of the model for shear flows, the new model is re- 
quired to reduce to small scale isotropy if its source 
vanishes. 

5. Correct one-point anisotropy. The combined effects 
of the new stochastic diffusion terms with 0j , 0^:, , 
hij and Gij should be the correct level of one-point 
anisotropy in the correct components of bij . 

6. Consistent turbulent kinetic energy. Introducing 
anisotropy at the small scales must not create spu- 
rious turbulent kinetic energy. In other words, the 
kinetic energy budgets of the VD PDF model and 
of the Navier-Stokes equation must be consistent 
with each other, similarly to Eq. (11) in constant- 
density flows. 



D. Summary 

The extension of the constant-density Langevin model 
to variable-density pressure-gradicnt-driven flows com- 



prises of (1) working with the instantaneous mass- 
density, (2) relaxation towards the Favre mean (instead 
of the Reynolds mean) velocity, and (3) accounting for 
small scale anisotropy with a new tensorial diffusion term 
(instead of built-in local isotropy) . 

In the next sections we examine the equations govern- 
ing the moments of the VD SDE (13). The derived mo- 
ment equations are rigorous mathematical consequences 
of the particular functional form of Eq. (13), which es- 
tablish constraints on the coefficients, Gij, hij, (j), and 
(pm that serve as a guideline for model specification. 



IV. REYNOLDS MOMENT EQUATIONS 

To investigate the new features of the SDE (13), and 
to aid the specification of its coefficients, the evolution 
equations for the first two moments of velocity, derived 
from the PDF model and the VD Navier-Stokes equation, 
are compared and examined. Here we discuss ensemble 
(or Reynolds) averaged statistics as they explicitly show 
the effects of the density fiuctuations. Sec. V discusses 
the Favre moments. 



A. The PDF transport equation 

This section gives the transport equation from which 
the equations governing the Reynolds moments are de- 
rived. 

The equivalent Fokker-Planck equation governing the 
Eulerian joint PDF of density and velocity, f{g,v;x,t), 
is derived from the VD Langevin model, Eq. (13) and 
Eq. (8), see Ref. 4, 



dl 
dt 



' dxi 



d_ 

dvi 

d_ r 

dv 
1 



2 dvidv 



{,9i ~P,i/9 + t^/o^'^Vi) f 

[Vj))f 



G^j [vj 



(21) 



<^oH,,)ef 



+ density model terms, 

where Hij is a symmetric non-negative semi-definite ten- 
sor with its Cholesky-decomposition 



H, 



hibh 



(22) 



Beside specifying the coefficients, Gij, Hij, 4>i and 4>d, 
Eq. (21) can be closed if a model for the dissipation rate, 
e, is specified. Such a model is given by van Slooten et 
al.34 

The particular form of the density equation, that yields 
additional fluxes in the density sample space in Eq. (21), 
is not important for our purposes. A density model, 
based on the beta PDF, is given in Ref. 2. However, 
the beta model is only one of many possibilities and thus 
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we do not detail the functional form of the density pro- 
cess, Q*{t). Eq. (21) only signals the availability of the 
instantaneous density field (as a sample space variable, 
g) and thus the full density PDF. 

In PDF methods, terms originating from the physi- 
cal process of advection appear in closed form, thus we 
will incorporate these in the Lagrangian derivative. As 
a consequence, the derived moment equations will repre- 
sent the rate of change along instantaneous Lagrangian 
paths, as specified by 



d 



d 



At dt ^'^ dxk 



(23) 



This allows us to concentrate on the fundamental physics 
due to pressure-gradient and viscous forces and their 
modeling, separately from advection. It is worth empha- 
sizing that the terms underlying d(-)/d<: may be different 
in equations governing different statistics. 

In a joint PDF formulation for density and velocity, 
both Favre and Reynolds-averaged statistics can be ob- 
tained from the PDF. Variable-density flows are tradi- 
tionally investigated in the Favre-averaged framework, 
in which the moment equations take a simpler form at 
the price of concealing some effects of the density fluctu- 
ations. We work here with Reynolds-averaged statistics 
for three reasons: 

1 . With Reynolds averaging the effects of density fluc- 
tuations are made clear. 

2. All new features of the extended model can be in- 
vestigated through equations for the Reynolds mo- 
ments, as will be seen. 

3. In the PDF framework modeling is not performed 
explicitly on the moment equations, nor are these 
explicitly discretized and solved. Therefore, from 
the viewpoint of analyzing the features of the veloc- 
ity model, it is unimportant which type of statistics 
we investigate. 

In the following, we concentrate on the new features of 
the extended model, Eq. (13), in VD flows where the 
density fluctuations are large, compared to the GLM, 
Eq. (7), in constant-density flows. The Favre-averaged 
moments and the Boussinesq limit are discussed in the 
subsequent sections. 



B. The ensemble mean velocity equation: Vi 

The modeled Reynolds-averaged mean velocity equa- 
tion in VD flows is derived from Eq. (21), see Ref. 4: 



dt 



gi - v-p,i +fiv-\/'^Vi + G^-jv'^, (24) 



SDE 



with the Favre fluctuation of the velocity, d" — Vi — {vi 



For comparison to the constant-density case, the corre- 
sponding mean velocity equation, derived from the gen- 
eralized Langevin model, Eq. (7), reads 



dvi 
'dt 



9i - Pn /qo + fJ-fga^^Vi 



(25) 



GLM 



Comparing Eq. (24) to its counterpart, derived from 
the Navier- Stokes equation for VD flows, Eq. (3), 



"dF 



gi — v-p,i +fiv-V Vi — v'p',i + /iw'V^f^, (26) 



NS 



indicates that the combined effects of the fluctuating 
speciflc volume (last two terms) are represented by the 
variable-density model as 



v'p',i + /iu'V^Uj' 



(27) 



The comparison of the right hand sides of Eqs. (24) 
and (26) also shows that employing Vj instead of (vj) 
in Eq. (13) would diminish the effect of v' on the mean 
velocity, since it would imply —v'p',i-\- fj^v'S/^v'^ = 0. This 
is not justified in VD flows with large density fluctua- 
tions. In contrast, Eq. (27) indicates how the effects of 
the fluctuating specific volume are incorporated into the 
velocity model for VD flows. 

To summarize, the product of the tensor G,j and the 
mass flux, w", jointly represents the correlations of the 
speciflc volume with the pressure gradient and viscous 
diffusion, Eq. (27). Compared to the constant-density 
case, this is an additional task, required of Gij, in VD 
flows with large density fluctuations. 



C. The ensemble Reynolds stress equation: v[v'j 

The model equation for the ensemble Reynolds stress, 
v'^Vj, in VD flows, derived from Eq. (21), becomes 



dt 



SDE 



+ {4>i6ij + 4'nHij)e, 



(28) 



where the identity v[v'^ = v[v'j has been used in the terms 
involving dj. For comparison to the constant-density 
case, the GLM equation (7) yields 



dv'jv'. 



dt 



(29) 



GLM 



In Eq. (28) the tensor Mij is defined as the integral 
over the whole sample space of the joint density- velocity 
PDF," 



(-lJJ3,.j +/.wV^t)i) / 



dgdv 



= —v'v'--p,j —v'v'--p,i +iJL{y'v[-V'^Vj -I- w'w^- • V^Wi) . 

(30) 
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Aiij collects the effects of the specific volume flux, v'v'^, 
due to the mean pressure gradient and the mean viscous 
forces. The products —v'v'^-p,j —v'v'^ -pji are an important 
source of turbulence in VD flows and reflect the fact that 
Lagrangian particles of different-density fluids accelerate 
very differently in response to pressure gradients. 

Since gv = 1, the speciflc volume flux, v'v'j, appearing 
in the Reynolds-averaged formulation, is related to the 
mass flux by 



Q-v'v[ + v-q'v[ + q'v'v[ 



0. 



(31) 



Consequently, v'v[ can be expressed in terms of q'v[, and 
Mij can be written as 



{g'v'v'.-p,.j +Q'v'v'.-p,i)/ g 
fJ,{g'v'vl-\7'^Vj + g'v'v'y\7'^Vi)/g, 



(32) 



where at = g'v[/'g. All one-point statistics of p, v and 
Vi can be extracted from the joint PDF, /, therefore 
the terms in Eq. (32) require no explicit closure assump- 
tions in the PDF formulation. The representation of M.ij 
is consistent, provided the joint density-velocity PDF is 
valid. The necessary and sufficient conditions that estab- 
lish the validity of the PDF, /(p, v\ x, t), are:'' 



1. Realizability: 



f{g,v;x,t) > 0. 



2. Normalization property: 



f{g, v; X, t)dvdg = 1. 



(33) 



(34) 



3. Conservation of mean mass: f(g,v;x,t) satisflcs 
conservation of mass in the mean. 



^ = -g-d- g'd\ 



(35) 



with d 



In other words, if a mass-density model, satisfying the 
above three conditions, is coupled to the VD SDK (13), 
the representation of Mij remains consistent. Such a 
density model is discussed in Ref. 2. 

Eq. (30) indicates that Aiij vanishes in constant- 
density flows, as v' = 0. Consistently, Aiij = in the 
constant-density model, Eq. (29). Its appearance in the 
VD Reynolds stress budget, Eq. (28), is the consequence 
of employing the instantaneous particle density, g*, in 
the mean forces terms of the SDE (13). 

The Favre Reynolds stress equation, derived from the 
SDE (13) and the Navier-Stokes equation (3), are dis- 
cussed in Sec. V. As will be shown, in that framework 
the effects of the mass flux in the first line of Eq. (32) 



appear explicitly, while the terms involving triple correla- 
tions are embedded in the Favre Reynolds stress, {v'/v'J). 

The equation governing the ensemble Reynolds stress 
in VD flows, derived from the Navier-Stokes equation (3), 
is 



dt 



(36) 



NS 



with the redistributive, transport, dissipative terms, and 
the triple correlations, respectively. 



Tij 



V-p'{v'ij + V'j,i), 



{v[p'),,+{vy)n~pv''v^ 



(37) 
(38) 

(39) 

V^J = -v'{v[p',j +v'jp',, ) + fiv'iv'y^v'^ +v'jV^v'^). (40) 

Note that production due to mean deformation, "T^ij — 

v[v'f.-Vj^k+VjV'^.-Vi^k, and turbulent transport, iv'yjV'f.),k, 
remain hidden in the above development. These terms 
arc in the Lagrangian derivative, dv'^v'j/dt. 

Comparing the model (28) and Navier-Stokes (36) 
Reynolds stress equations in VD flows shows that the 
terms involving Gij and e provide a joint model as 



V,, 



(41) 



Eq. (41) indicates that the terms in Gij and e in the 
SDE (13) represent the combined effects of the pressure- 
strain correlation, turbulent transport, dissipation and 
the triple correlations. It also shows that A^^j does not 
require modeling. The model terms of Eq. (13), the 
ones involving Gy and e, do not account for M.ij, which 
appears naturally in the moment equation (28). M-ij 
is the most important term that drives variable-density 
pressure-gradient-driven turbulence. 

In second-order moment closures for VD flows equa- 
tions for the first two Favre moments, (wj) and {v'/v'J), 
are solved. Since only Favre-averages are involved, it 
is necessary to separately model the normalized mass 



fiux, Oi 



to compute its effect on the Reynolds 



stress. One such model is by Besnard et al.,"^^ which car- 
ries two model equations in addition to the Favre mean 
and Reynolds stress: one for Oi and another one for the 
density-specific- volume covariance, b = ~g'v', appearing 
in the equation for the mass flux, each with its burden of 
closure assumptions. As Eqs. (30) and (32) show, in the 
PDF formulation no modeling of a; and b is required. If 
the joint PDF is valid, the representation of both and 
b are consistent, see also Ref. 2. In the PDF framework 
modeling is performed on the joint PDF (instead of its 
moments), e.g. on the stochastic equations representing 
g* and v* . In the case of the SDE (13) only the effects of 
the fluctuating pressure gradient and viscous dissipation 
arc modeled, while the effects of the mean forces, from 
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which originates, are exact, see also Eqs. (17) and 
(41). 

To summarize, in contrast to moment closures, the im- 
portant source of turbulence in VD flows, Jvi^j Eq. (32), 
is closed in the PDF formulation. Redistribution, tur- 
bulent transport, dissipation and the triple correlations 
of the specific volume and the gradients of velocity and 
pressure are jointly modeled, Eq. (41). 



D. The ensemble turbulent kinetic energy 
equation: k 

Taking one-half the trace of the Reynolds stress equa- 
tion gives an equation for the turbulent kinetic energy. 
Contracting Eq. (28) results in the model equation for 
the Reynolds-averaged k = v[v[/2 in VD flows as 



the GLM Eq. (7) creates no spurious turbulent kinetic 
energy in homogeneous turbulence. In contrast, Eq. (49) 
is a model for variable-density inhomogeneous flows with 
small scale anisotropy. As the joint PDF, f{g,v), con- 
tains no information on the fluctuating pressure, its gra- 
dient or the velocity derivatives, their correlations require 
modeling. Eq. (49) indicates how the specifications of 
Gij, Hij, (jjj and 0^ yield a combined model for the ef- 
fects of the pressure-dilatation covariance, the trace of 
transport and the triple correlations, TZ — T + V. TZ is 
important at high turbulent Mach numbers, while V at 
high Atwood numbers, see also Sec. VI. 

To summarize, Eq. (49) provides a model constraint 
in inhomogeneous variable-density flows with small scale 
anisotropy for the effects of transport and scalar effects 
of variable density. 



dfc 
dt 



M + G,M + 



SDE 



(42) 



E. The ensemble Reynolds stress anisotropy 
equation: bij 



with A4 = M.ii/2 and H = Hii/2. For comparison to 
constant-density flows, the GLM Eq. (7) produces the k 
equation 



dfc 
dt 



GLM 



(43) 



From Eq. (36), the equation governing fc based on the 
Navier-Stokes equation (3) is 



dfc 
dt 



M+TZ-T -e + V, 



NS 



with the half of the traces of Eqs. (37-40), 
T=v 



e ■ 
V 



-v'v[p',i -\- fiv'v^V^v' 



(44) 



(45) 
(46) 

(47) 
(48) 



Comparing Eqs. (42) and (44) shows that both yield the 
same evolution of fc if 



G^vlvr +1^1 + -<t>: + c^^Hj e = n-T+V. (49) 

Eq. (49) can also be obtained from contracting Eq. (41). 

As Aiij appears closed, its trace, A4, requires no ex- 
plicit modeling. Eq. (49) in inhomogeneous VD flows 
with small scale anisotropy can be contrasted to Eq. (10) 
in homogeneous constant-density flows where the small 
scales are assumed to be isotropic. The differences are: 
the trace of the anisotropic diffusion coefficient tensor, 
H, pressure-dilatation, 7?,, pressure and viscous trans- 
port, T, and the trace of the triple correlations, V. 

In constant-density homogeneous flows 7?, = T = V = 
and Eq. (10) is a consistency condition that ensures that 



The normalized one-point Reynolds stress anisotropy 



I 3 

2k 



(50) 



is indicative of the fraction of the turbulent kinetic energy 
in a given component of the Reynolds stress tensor: —1/3 
indicates no energy, while 2/3 indicates 100% energy in 
the component. Its importance in both constant-density 
and VD flows can be highlighted by writing the produc- 
tion term of v'-v'i due to mean deformation as 



V^J = 2fc ( b,k 



^^ikj V],k + '^k (bjk + ^Sjk ) Vi^k- (51) 



In our formulation, Vij is hidden in the Lagrangian 
derivative of the Reynolds stress balance Eq. (36), in as 
much as 



dt 



NS 



dt 



Vk 



k 



dv[v',jV 



dxk 



(52) 



Eqs. (51) and (52) show that the correct prediction of 
both isotropic and deviatoric parts of v'^^ are important 
in representing the correct turbulence levels. To inves- 
tigate how anisotropy is created and dissipated by the 
stochastic VD model, the evolution equation governing 
hij is investigated in the following. 

The general case. The hij equation of the VD PDF 
model is derived based on Eqs. (28) and (42): 



2fc^ 
dt 



Mi + 2fc 



SDE 



1 



Gikhjk + Gjkhik + "^(Gij I ^ji) 



2{M + Guv'yi) + (30, + 20^i/) 



Gkiv'^vlS^j +(l>neH^j, 



(53) 
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with the deviatoric parts of Mij and Hij denoted by the 
superscript "d" : 



H, 



(54) 
(55) 



Eq. (53) is the most general ease of 6y that ean be repre- 
sented by the VD SDE (13): it assumes variable density 
Mij 7^ 0, anisotropic Gij and anisotropic small scales 
Hij 7^ 0. Depending on the approximations that can 
be made, several simplified forms of Eq. (53) can be ob- 
tained. 

(1) GLM in constant-density flows. As a special 
case, the GLM Eq. (7) for constant-density flows yields 



2k ^ 
dt 



2k 



GLM 



Gikbjk + Gjkbik + '^i^ij + ^ji) 



- (2Gkiv'kvi + 3Co£) - ^Gklv'^v'^6^J. (56) 

This amounts to applying the GLM in constant-density 
flows, where Gij is given by an anisotropic model, as 
discussed in Pope.^^ It can be obtained from the general 
case by assuming constant density Mij = 0, anisotropic 
Gij and isotropic small scales Hij = 0. 
(2) SLM in constant-density flov^rs. Eq. (56) further 
simplifies if Gij is chosen to be isotropic, e.g. in case of 
the SLM, Eq. (12), since all terms involving Gij vanish, 
as expected: 



dbjj 



SLM 



2^0 -^bij, 



(57) 



resulting in Rotta's well-known linear return-to-isotropy 
model,i'27 if Co = 2/i{GR - 1), where Cr is Rotta's 
constant. This can be obtained from the general case by 
assuming constant density AAij = 0, isotropic Gij and 
isotropic small scales Hij ~ 0. 

(3) GLM in VD flows with small scale isotropy. If 

Gij is chosen to be anisotropic, while the small scales to 
be isotropic, i?.y=0, Eq. (53) simplifies to 



2k ^ 
dt 



2k 



GLM 



Gikbjk + Gjkbik + 2 (^ii + ^ji) 



2(^M+ Gfc/<v;) + 3(f>re 



Gkiv'^v'iSij, 
(58) 



where the superscript "i" indicates presumed small scale 
isotropy. This case amounts to applying the GLM in 
VD flows, where Gij is given by any anisotropic model, 
such as discussed by Pope.^® It can be obtained from 
the general case by assuming variable density Mij ^ 0, 
anisotropic Gij and isotropic small scales Hij = 0. 
(4) SLM in VD flows with small scale anisotropy. 
If an isotropic Gy is chosen to be employed in the ex- 
tended model, Eq. (13), for VD flows, in Eq. (53) the 



terms involving Gij cancel as in the constant-density case 
and bij will be governed by 



2k^ 
dt 



= Mf,+ e(f>nH° 



SLM 



(59) 



- [2M + {3(br + (boHkk)e] b,j, 

where the superscript "a" denotes presumed small scale 
anisotropy. This can be obtained from the general case 
by assuming variable density Mij 7^ 0, isotropic Gij and 
anisotropic small scales Hij ^ 0. 

(5) SLM in VD flows with small scale isotropy. 

Eq. (59) further simplifies in a case that adheres to Kol- 
mogorov's hypothesis of local isotropy, i.e. Hi 



0, 



2fc^ 
dt 



i2M + 30,£) 6, 



(60) 



SLM 



This can be obtained from the general case by assuming 
variable density Mij ^ 0, isotropic Gij and isotropic 
small scales Hij ~ 0. 

In constant-density high-Rcynolds-number flows that 
are free of large strains and rotations the small scales are 
isotropic. In contrast, Rayleigh- Taylor flows are aniso- 
tropic at all times at both large and small scales. ^"^^'^'^'^^ 
These variable-density flows are at moderate Reynolds 
numbers and free of large distortions. To this end we are 
concerned with both bij and dij. 

In the homogeneous RT DNS of Livescu & 
Ristorcelli,^'^° Vij = thus Mij is the sole source of 
bij which quickly vanishes in the initial stages of the flow 
evolution. On the other hand, 6^ indicates anisotropy for 
the full extent of the simulations. This is most likely due 
to the suddenly decreasing turbulent Reynolds number 
closely following the behavior of Mij. As the signifi- 
cance of the non-linear term in the Navier-Stokes equa- 
tion (compared to viscous dissipation) decreases, the ef- 
ficiency of the non-linear cascade, that would transfer 
energy from the large to the small scales, weakens. This 
effectively locks in the anisotropy structure of bij which 
thus dissipates very slowly. 

The VD PDF model is designed to account for a source 
of anisotropy that is independent of Mij and a (possibly) 
initially nonzero bij. That is, in principle, the model is 
capable of capturing the large and small scale anisotropy 
of variable-density, strained or low-i?e flows through the 
combined effects of Gij and Hij. For generality, the fol- 
lowing discussion supposes an initially zero bij. 

The special cases of Eq. (53) above show that there are 
two independent sources of anisotropy in the VD velocity 
field, represented by the SDE (13) (without significant 
shear-production): via Gij and via Hij. This is indicated 
by 



1. Eq. (58), where H^^ 



and the source is 



proportional to the anisotropic Gij, (large scale 
anisotropy), and, 

2. Eq. (59), where the isotropic Gij has no effect 
and the source is proportional to Hij, (small scale 
anisotropy). 
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In the absence of shear production, Vij = 0, embedded in 
the Lagrangian derivative, Eq. (60) prescribes a decay of 
hij after vanishes. If we assume Hij to be nonzero 
we can take two routes in defining Gij : 



1. Isotropic Gij, special case (4)- In this case bij is 
governed by Eq. (59). If one chooses an isotropic 
Gij in VD flows, Eq. (59) shows that a nonzero 
Hij can sustain 6y- independently even after Mij 
vanished, c.f. Eq. (57). An isotropic Gij also means 
that the effects of large density fluctuations on the 
mean velocity are assumed to be isotropic, since 
according to Eq. (27), ~v'p',i + fiv'\7'^v'^ = -vf/T, 
where T denotes some function of a time scale, e.g. 
k/e. 

2. Anisotropic Gij , the general case. In this case bij is 
governed by Eq. (53). The choice of an anisotropic 
Gij combined with a nonzero Hij opens up all the 
possibilities of the SDE (13) to correctly capture bij 
in VD flows. This most general formulation seems 
to be capable of representing nonzero bij patterns 
that are more complex than that of a homogeneous 
RT flow. 

In summary, we showed that the VD SDE (13) allows 
for two independent sources of anisotropy, proportional 
to Gij and Hij, respectively. The question of whether 
Gij , which is capable of sustaining anisotropy in the ab- 
sence of shear production, should be assumed isotropic 
or anisotropic in general VD flows remains an open ques- 
tion. Functional forms of Gij and Hij are the subject of 
the next paper, Ref. 21. Their combined effect should 

1. Create the correct 6y structure in VD flows by 
modeling the processes in the Reynolds stress bud- 
get as Eq. (41) prescribes, 

2. Gij should provide the correct representation of the 
processes 



-v'p',i 



velocity equation (24), and. 



Gijv'J in the mean 



3. Gij and Hij should provide the correct representa- 
tion of the variable-density effects on the turbulent 
kinetic energy as given by Eq. (49). 



F. The functional forms of (f), and (po 



of Eq. (61) affects the turbulent kinetic energy budget, 
while only the deviatoric part affects the bij budget. This 
is ensured if 



b, = Co--Hcl)„, 



(62) 



which can be easily seen from the corresponding terms 
in the k equation (42), 



and the bij equation (53), 



^H ] e ~ 2^0^, 



(3(/), + 2cf>^H)eb,, = SCoeb, 



(63) 



(64) 



Eqs. (63) and (64) respectively show that with the choice 
of (f), in Eq. (62), k will not be affected by Hij or 
<j)^, while bij will only be affected by Hf^ and 0d, see 
Eq. (53). Eq. (62) ensures that the VD model reduces to 
its constant-density counterpart with small scale isotropy 
if — >■ 0: Eq. (42) governing k approaches Eq. (43) in- 
dependent of Hij . 

Since 4>i must be non-negative, the following bounds 
on (po must also be imposed: 



- ^ 2H 



(65) 



As Co and Hij are bounded, Eqs. (62) and (65) ensure 
that and 0d are also bounded. Thus we specify (j)^ as 

0o = [1-5(0)]^ with Q<g{e)<l, (66) 

where giO) is some function of the mixing state, 0, with 
g{9) = in the unmixed state and g{0) = 1 in the fully 
mixed state. The simplest specifications for g{0) and 
are through the first two density moments*'^^'^^''^^ 



9{e) ^e = i-- 



q{1 - q)' 



(67) 



Since g'"^ is maximum in the unmixed state and vanishes 
in the fully mixed state, Eq. (67) ensures the required 
behavior of g{9) in the extreme states. 

Incorporating the new specifications in Eqs. (62) and 
(66) for (t>i and (/)d into Eq. (13), the final form of the 
VD SDE is 



This section establishes further constraints on the func- 
tions (f)i and (j)o , governing the diffusion terms in the VD 
SDE (13), 



.1/2 



dW,: 



(61) 



d< = (gz - Pn /g* + fi/g*V^Vi) dt + {v* - (w^)) dt 



[G,gie)ef'dW.. + \'^ 



1/2 

[1 - g{e)]e } h^jdw;. 



As a reminder, the mathematical representation of a 
diffusion process requires cj^j, (jy^) and hij to be non- 
negative and bounded, which also ensures realizability. 
From the modeling point of view, it is customary to 
construct these functions so that only the isotropic part 



where Hij = hikhkj and H = Hii/2. Eq. (68) shows 
that Hij can have any unit, as it only enters in a nor- 
malized fashion and that the only restrictions on its 
functional form arc being symmetric non-negative semi- 
definite and bounded so that statistics of the process exist 
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and arc well defined. Substituting Eqs. (62) and (66) into 
Eqs. (28) and (42) gives the final governing equations for 
the Reynolds stress and turbulent kinetic energy, respec- 
tively, 



dt 



SDE 



dfc 
dt 



SDE 



(69) 
(70) 



These equations show that only the normalized devia- 
toric part of Hij affects the second moments and only 
their anisotropic part, i.e. hij. As mixing progresses g{9) 
takes on values from to 1 and diminishes the small scale 
anisotropy, independent of Hij . 

The developed constraints on the diffusion coefficients 
4>i and 4>D, Eqs. (62) and (66), ensure that (1) only the 
isotropic part of the sum of the two stochastic diffusion 
terms affect the evolution of k and (2) hij is affected by 
only the normalized deviatoric part. Enforcing mathe- 
matical consistency on the constraints, a function of a 
mix metric, g{9), naturally appears which ensures that 
the small scale anisotropy is diminished in the fully mixed 
state. The SDE (68) reduces to its constant-density coun- 
terpart, Eq. (7), as q' — > 0, i.e. g* — > and g{9) — > 1. 



This allows the equation to capture both large {hij) 
and small {dij) scale anisotropy of variable-density, 
strained or low-i?e flows. 

• We developed constraints on the diffusion coeffi- 
cients, Eqs. (62) and (66). These ensure that the 
isotropic part of the stochastic diffusion terms in 
Eq. (13) only affects the turbulent kinetic energy, 
while the Reynolds stress anisotropy is only af- 
fected by the deviatoric part. A mix metric ensures 
vanishing small scale anisotropy of different-density 
species in the fully mixed state. 

• The final form of the VD SDE (68) (incorporating 
the constraints on the diffusion terms) reduces to 
the GLM for constant-density flows as g' — >■ 0. 



V. FAVRE MOMENT EQUATIONS 

Variable-density flows are traditionally investigated 
(and modeled) using Favre-averaged variables. This sec- 
tion derives the Favre moment equations of the VD veloc- 
ity model, Eq. (13), and contrasts them to their counter- 
parts based on the Navier-Stokes equation. This is useful 
in comparing the PDF formulation to existing Favre mo- 
ment closures and provides insight for modeling unclosed 
terms. 



G. Summary of the ensemble moment equations 



A. The MDF transport equation 



We derived and analyzed the first two moment equa- 
tions governed by the SDE (13). The findings can be 
summarized as follows: 

• The drift, Gij{v* — (uj)), affects the Reynolds mean 
velocity in VD flows. This must be taken into ac- 
count by the specification of as indicated by 
Eq. (27). 

• In contrast to moment closures, the important 
source of turbulence in VD flows, Mij defined by 
Eq. (30), appears closed in the PDF formulation 
and thus requires no explicit modeling. The rela- 
tion (41) guides the modeling of the physical pro- 
cesses of redistribution, transport, dissipation and 
the triple correlations. 

• We derived a model constraint in inhomogeneous 
variable-density flows with small scale anisotropy, 
Eq. (49), for transport and scalar effects of variable 
density. 

• The VD SDE (13) allows for two independent 
sources of anisotropy, proportional to Gij and Hij , 
respectively. Gij is responsible for anisotropy at 
the large scales, while Hij at the small scales. 



This section gives the transport equation from which 
the equations governing the Favre moments arc derived. 

The transport equation governing the joint mass den- 
sity function (MDF), 



.^{g,v;x,t) EE gf{g,v;x,t), 



(71) 



can be obtained by multiplying Eq. (21) by g and using 
the law of mass conservation. 



dxi 



_d_ 

dv, 
d_ r 
dv 
1 



-p,^/g + ^l/g\/'^v^)^ 



2 dvidvi 



{cj)j5ij + cj)^Hij)e^ 



density model terms. 



(72) 



where for simplicity, we kept (pi and (po in the diffusion 
terms. This equation is integrated to obtain Favre mo- 
ment equations. In Favre averaging the density acts as a 
weight, therefore the precise functional form of the den- 
sity equation is unimportant from the viewpoint of de- 
riving Favre moment equations of the velocity fleld. Two 
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important characteristics of ,^ are: 

■^{g,v]x,t)dQ(lv —'Q{x,t), (73) 
Q{v).^{Q,v;x,t)<lQ<iv = g{x,t){Q{x,t)), (74) 



where Q{v) is almost any function'* with its Favre average 
(Q) . These attributes of ^ serve as the basis of obtaining 
the Favre moment equations. 



B. The Favre mean velocity equation: (ii^) 

Muhiplying Eq. (72) by the sample space variable Vk 
and integrating each term^ yields the equation governing 
the modeled Favre mean velocity, {vi), governed by the 
SDE (13), as 



dg{i 



at 



SDE 



dg{vt){vj 
dxi 



= Q9i-Pri +A*V 



(75) 

In moment closures {vi) is typically written as Ui or Vi. 
Contrasting Eq. (75) with its counterpart, derived from 
the variable-density Navier-Stokes equation (3), 



dQ{vi 



at 



NS 



dg{vi){v-j 



dg(v'-v' 



"v") 



Q9i 



(76) 

shows that the SDE (13) is consistent with the Navier- 
Stokes Favre mean. All terms are the same, in closed 
form, containing no explicit model terms. 

For comparison and later use, the Reynolds-averaged 
mean velocity equation, derived from the Navier-Stokes 
Eq. (3), is given by Eq. (26), or equivalently. 



'dt 



NS 



(77) 



where the effects of the fluctuating specific volume (last 
two terms) are exposed and have been shown to be mod- 
eled as —v'p'yi H- fiv'\7'^v[ = GijVj. 

The above development shows that the PDF model 
Eq. (13) is consistent with the Navier-Stokes equation 
in both Favre and Reynolds-averaged frameworks in the 
mean. The effects of large density fluctuations appear 
hidden in the former and explicitly modeled in the latter. 
This is in contrast to Favre moment closures, where they 
are not accounted for if a constant-density model is used 
to represent the Favre Reynolds stress, (v-'w"). 



and Eqs. (76) and (26) from Navier-Stokes: 

_dai 

SDE 



dt 



-g'v' {p,i -fiV^v,) + gGijUj 



(78) 



_dai 



= - g'v' {p,i -/xV^Ui) + g{v'p',i - fiv'V^v'^. 



NS 



Comparing Eqs. (78) and (79) yields 



- v'p',i + ^w'V^v^ = -GijUj, 



(79) 



(80) 



the same as Eq. (27), already obtained from contrasting 
the Reynolds mean velocity equations. Eq. (80) estab- 
lishes the consistent modeling of the mass flux by the 
PDF model and highlights the importance of GijUj as a 
model for the specific-volume-pressure-gradient and the 
specific- volume- viscous- force covariances. 

The above development can be put in context with 
VD Favre moment closures in the Eulerian framework 
by expanding the Lagrangian derivative in Eq. (79) to 
yield the mass flux equation*^ 



.da,; 
dt 



dTgai dg{vj)ai 



NS 



dt 



dxi 



9{a.a,),,-~^{V^--g{v'/v';))+^^ 



= -q'V {pn -^lV^v^) + g{v'p',, - iiv'V^v'^). (81) 

In moment closures, equations arc solved for p, {vi), 
{v'lv'^), b = -Yv' (or ^) and a,. In Eq. (81) the ef- 
fects of the unclosed terms, v[d' , 'g,jg'v'^v'p ig'v'^v'j),j and 

v'p',i — /iw'V^w^ are modeled. In contrast, in the PDF 
method, Eq. (80) shows that only the non-Boussinesq 
effects, v'p',i — pv'V^v'^, require closure assumptions. 

It is useful to investigate the model mass flux and mean 
velocity equations in the limit Vi = 0. This holds exactly 
at arbitrary Atwood numbers in a homogeneous RT layer 
and at low Atwood numbers in the inhomogeneous RT 
layer. Also, Vi w is still a good approximation in high- 
Atwood-number inhomogeneous RT layers.^ If Vi = 0, 
then Qi — (vi) , thus from Eqs. (75) and (78) we have 



- v■p,^ -(w'w-),j +vid' = Gi^ 



(82) 



The same is obtained by setting = in Eq. (24) 
with the advection term expanded as in Eq. (77). While 
Eq. (80) holds in inhomogeneous VD flows with large 
density fluctuations, Eq. (82) can be considered as an 
approximation. The importance of Eq. (82) is that its 
processes may be easier to model than those in Eq. (80). 



C. The mass flux equation: ai — q'v'^/ g 



D. The Favre Reynolds stress equation: {v'i'v'j') 



The equations governing the mass flux can be obtained 
from Qi = {vi)—Vi, using Eqs. (75) and (24) for the model 



Multiplying Eq. (72) by {vk ~ (wfe>)(w; - (u;)) tl^cn inte- 
grating produces the model equation governing the Favre 
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Reynolds stress, of the SDE (13) as 



dt 



SDE 



dxk 
—/„,//, //\ 



dxi, 



+ gGik {v'-v'D + gGjk {viv'D + 'g{(t)i5r 



(83) 



The Favre Reynolds stress equation, derived from the VD 
Navier-Stokes Eq. (3), is 



dt 



NS 

+ g{vtVk) {vj),k +g{v''vk) {v,\k = 



dxk 



dxh 



(84) 



+ {n^J -%j -eij)/v. 



Comparing Eqs. (83) and (84) shows that the first three 
lines are the same. The terms proportional to rep- 
resent the source of turbulence in the Favre-averaged 
framework whose role played by Aiij in the ensemble- 
averaged framework, c.f. Eqs. (36) and (84) and see 
Eq. (32) that details the effects of the fluctuating spe- 
cific volume. The correspondence of the remaining terms 
produces the modeling constraint 



(85) 



{Tit, - %j - eij)/{g ■ v) = 
= Gik{v'-v'l) + Gjkiv'lv'l) + {4>i6ij + (j)oHij)e, 

which can be contrasted to the one obtained from the 
comparison of the ensemble Reynolds stress equations, 
Eq. (41), to yield an expression for the triple correlations 
in Vij as 



V^J = -V'{v[p\j +V',p'„ ) + pLV'(v['^'^v'. + V2u^) 



g'v' 



g-v{Gikajak + GjkaiUk) ■ 



(86) 



Similar to the ensemble-averaged framework, the half of 
the trace of Eq. (85) produces 



{n-T-e)/{-g-v) = G.,{v':v'') + 



(87) 



A set of modeling constraints, based on the Favre 
Reynolds stress and turbulent kinetic energy, Eqs. (85) 
and (87), has been derived. These are analogous to 
Eqs. (41) and (49) in the ensemble framework and do 
not explicitly involve the triple correlations, Vij. 



E. 



Summary of the Favre moment equations 



We analyzed the Favre velocity moment equations, de- 
rived from the PDF model, that require modeling in VD 
moment closures. The main findings of this section can 
be summarized as follows: 



• The PDF and Navier-Stokes Favre mean velocity 
equations are consistent. The effects of the fluc- 
tuating specific volume, hidden in the Favre for- 
mulation, are explicitly modeled in the ensemble 
framework, discussed in Sec. IV. 

• We derived a modeling constraint that relates the 
model tensor Gij to the effects of the fluctuating 
specific volume, Eq. (82), which is as an approxi- 
mation for Eq. (80) for flows with Vi « 0. 

• We showed that the source of turbulence in VD 
flows in the Favre-averaged framework (the effect of 
the mass flux, a^) appears closed in the PDF formu- 
lation. Pressure redistribution, transport and dissi- 
pation of the Favre Reynolds stress tensor and tur- 
bulent kinetic energy are jointly modeled according 
to Eqs. (85) and (87), respectively. These two equa- 
tions can be used as a guideline for specifying Gij 
and Hij if the Favre framework is preferred. 



VI. THE BOUSSINESQ LIMIT 

The well-known Boussinesq approximation, in which 
the effects of variable density arc only retained in the 
body force in the Navier-Stokes equation, is used in many 
situations. It is necessary to ensure that the VD SDE 
(13) can account for the Boussinesq limit. 

A. The instantaneous Navier-Stokes equation: Vi 

To introduce the notation the VD and Boussinesq 
Navier-Stokes equations are given. 

Wc start from two separate forms, written for the 
variable-density case and the Boussinesq case, respec- 
tively. 



di'i 
'dT 

dt 



VD 



go' 



Q9i - Pli Vi, 



ggi -p,i+^i^ Vi, 



(89) 



where d{-)/dt denotes the derivative along an instanta- 
neous Lagrangian path. The density on the left hand side 
of Eq. (88), with the aid of the specific volume, g = 1/v, 
can be decomposed as g = 1/iv + £_v'). Thus the above 
two equations can be written in a combined form that ex- 
plicitly shows the departure from the Boussinesq-limit: 



dvj 
'dt 



g{v + ^v')g,- {v + ^v')p„+p{v + ^v')V^v,. (90) 



The non-dimensional ordering parameter, ^= [0, 1], dis- 
tinguishes between the Boussinesq case with ^ = 0, when 
g = go ='g = 1 /v, corresponding to 



dvi 
'dt 



gvg^ - vp,i +/iw-V Vi, 



(91) 
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and the variable-density case with ^=1, when g-v+g'v' 
1, corresponding to 



dF 



gi - vp,i +fiv'V Vi, 



(92) 



VD 



the same as Eq. (3). 

The equations governing the first two ensemble mo- 
ments of Eq. (90) are investigated in the next sections. 
Then they will be contrasted to the VD PDF model in 
the Boussinesq and VD limits. 



with TZij, Tij, Eij and Vy defined by Eqs. (37-40) and 



-l]; It] 



M 



Q{v'Vjgt + v'v[gj) 



VD 



+ g'v'v'^gi + Q'v'v[gj + M^J, 



(98) 
(99) 

(100) 



As before, setting ^ = and g ~ 1/v gives the Boussinesq 
limit, 



B. The Navier-Stokes ensemble mean velocity 
equation: Vi 

The equation governing the Reynolds-averaged mean 
velocity, derived from Eq. (90) is 



dv 
"dt 



g-vgi + v-p,i +fivV^Vt 



(93) 



+ ^{Q'v'gi - v'p'yi + ^v'V^v'^. 
Setting ^ = and 'g = 1/v gives the Boussinesq limit. 



dt 



gi - v-p,i +11V-V Vt, 



(94) 



showing that in Boussinesq flows, where v = const., the 
mean velocity along an instantaneous Lagrangian path is 
governed by the same equation as in incompressible flows 
if the density is chosen as go = ^ = 1 /v = const.. 



dvi 
'dt 



gi ~ p,i / go + iJ./ go^'^Vi. 



(95) 



Setting ^ = 1 in Eq. (93) gives the mean velocity equation 
in the VD case as 



dvi 
'dt 



^ gi - v-p,i Vi - v'p',1 + /iw'V^w-, (96) 



VD 



the same as Eq. (26), since 'g -v + g'v' = 1. Comparing 
Eqs. (94) and (96) indicates that the fluctuating specific 
volume is directly associated with large density fluctua- 
tions in VD flows, due to their non-negligible effect on 
the inertia terms and the departure from the Boussinesq 
limit. 



C. The Navier-Stokes ensemble Reynolds stress 
equation: v[v'^ 



dt 



'71 7 



(101) 



which does not involve the mean pressure gradient but 
contains the buoyancy force, gi, in Mfj, Eq. (98). Setting 
^ = 1 in Eq. (97) gives the ensemble VD Reynolds stress 
equation 



dvW, 



dt 



(102) 



VD 



the same as Eq. (36), which does not involve the buoy- 
ancy force but contains the mean pressure gradient, Pj, 
as the terms proportional to gi cancel in Aifj+AiYj^ due 
to the identity 'g-v'v'^ + v-g'v[ + g'v'v[ ~ 0. 



D. The VD Langevin model in the Boussinesq limit 

The asymptotic behavior of the SDE (13) for VD 
flows, as it reduces to the constant-density SDE (7) when 
g' — )■ 0, has been discussed in Sec. IV. Here we investi- 
gate the question: "Under what circumstances the mo- 
ments of the VD PDF model approaches the moments of 
the Boussinesq limit?" The procedure yields additional 
model constraints. 

The equations governing the ensemble-averaged mean 
velocity and Reynolds stress along instantaneous La- 
grangian paths, respectively, derived from the SDE (13), 
read 



dvi 
'dt 

dv[v'. 



gi - v-p,i +^iv-V'^Vi + Gi-jv'-, (103) 



dt 



VD 



VD 



(104) 



The equation governing the ensemble-averaged Rey- 
nolds stress, derived from Eq. (90) is 



dv[v'^ 
dt 



Ml + 7^., - - + ^(Al^D + Vy) , (97) 



the same as Eqs. (24) and (28). These equations consti- 
tute as models for Eqs. (93) and (97) in the general VD 
case with ^ = 1. 

Correspondence of the mean: TJ^. Comparing 
Eqs. (94) and (103) indicates that the VD PDF mean 
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velocity equation correctly reduces to the Boussinesq 
Navier-Stokes mean if the fohowing holds: 



[GijUj 



0. 



(105) 



Correspondence of the Reynolds stress: 

The asymptotic behavior of the model Reynolds stress 
equation can be most easily seen from the SDE that ex- 
plicitly shows the departure from the Boussinesq limit: 



dv* = {g*v*g, - v*p,i +fiv*WH,) dt + {v* - (vj)) dt 

(106) 



,e)'/' dW, + {cl)ney^''h,,dW', 



1/2 



where v* = v + ^w'. This is the same as Eq. (13) since 
V* = l/g*. The model Reynolds stress equation, derived 
from the SDE (106) becomes 



dt 



VD 



(107) 



£. 



In the Boussinesq limit, where ^ = and g = 1/w, 
Eq. (107) as a model correctly reduces to the Boussinesq 
Navier-Stokes Reynolds stress Eq. (101), provided 



GikvW^ + Gjkv'y^ + {(t),5ij + (t)nHrj)e = 

■' J B 



(108) 



i.e. if the effect of Vij in the VD model vanishes. Eq. (86) 
shows that this can be ensured if 



[GikUjak + Gjfcfljafcjg = 0. 



(109) 



In the fully VD case with ^ = 1 and g ■ v'v[ + v ■ g'v'^ + 
g'v'v[ — 0, Eq. (107) correctly gives the modeled Navier- 
Stokes Reynolds stress Eq. (102). 

We obtained two additional model constraints, 
Eqs. (105) and (109), to ensure the correct asymptotic 
behavior of the first two moments of the VD PDF model 
in the Boussinesq limit. The PDF model equations for 



and 



« 3 



can now be made consistent with that of the 



Navier-Stokes Eq. (90) with ^ = 0. 



E. Summary of constraints; VD and Boussinesq 
forms of the PDF model 



The set of equations, indicative of how the various 
physical processes are represented by the PDF model at 
the first two levels of statistical moments, have been sum- 
marized in Table I. This is the starting point for devel- 
oping possible functional forms of the coefficients, Gy , 
Hij and g{9). 

In the VD case the stochastic model equation is 



dvl 



VD = {9i - +^^v*\7'^Vi) dt + Gij {y* - {vj)) dt 



6^e)'/' h,,dW'. 



(110) 



Alternatively, if the development of a Boussinesq 
model is sought, the following equation automatically 
satisfies the requirements of the Boussinesq limit, Eqs. 
(105) and (109): 

dw*B = {g*'^9i ~ "^-Pri V^tJi) dt + Gij [v* - Vj) dt 



.1/2 



hudW', 



(111) 



whose first two moment equations can be easily seen to 



reduce to Eqs. (94) and the modeled (101) with V^. 



0. 



VII. CONCLUSION 

We developed and discussed fundamental ingredi- 
ents necessary for extending Langevin-type models from 
constant-density shear flows to variable-density pressure- 
gradient-driven turbulence. A forthcoming article^^ will 
combine these ideas with the ones on active scalar 
mixing^ to specify a joint probability density function 
(PDF) model for mixing-driven variable-density (VD) 
hydrodynamics and will present validation for Rayleigh- 
Taylor flows. 

In this paper we proposed a stochastic differential 
equation (SDE) for modeling the instantaneous veloc- 
ity increments of Lagrangian particles in VD turbulence. 
The functional form of the equation and its consequences 
have been discussed in detail. Several consistency condi- 
tions and constraints, based on mathematical and phys- 
ical considerations, have been developed. These results, 
summarized in Table I, are all rigorous mathematical 
consequences of the particular functional form of the pro- 
posed model equation and some basic physical consider- 
ations in VD flows. 

In particular, we showed that the widely used general- 
ized Langevin modcl^ for the velocity PDF in constant- 
density flows can be extended to VD flows to include the 
following features: 

1. Representing variable- density effects. The ex- 
tended model incorporates the effects of large den- 
sity fluctuations due to non-uniform species con- 
centrations on the fluid particle momentum in 
variable-density turbulence. This includes: 

(a) VD effects on the inertia and mean forces 
terms in closed form. These are cubic (as well 
as quadratic) non-linearities in the Navier- 
Stokes equation and the mixing arising from 
the strong coupling between the density and 
velocity flelds. 

(b) VD effects on the fluctuating pressure gradient 
and viscous forces in modeled form. 

2. Closed mass flux and density- specific-volume co- 
variance. In variable-density flows, these processes 
relate to a primary mechanism that generates tur- 
bulent kinetic energy through the mean pressure 
gradient. The joint PDF model for mixing-driven 
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TABLE I: Summary of main results. The PDF model for variable-density flows and its model constraints, involving the 
coefficients dj, Hij = hikhkj and g{d), related to various physical processes (defined at the right) at the first two levels of 

statistical moments. Note that the equation-pairs marked by (= and =) and {= and =) express equivalent constraints in the 
Reynolds and Favre-averaged frameworks, respectively. 

c\v: = {g^ - P,^ Iq* + M/e* V'«,)df + G., [y] - {vf)) At + {4'iet''' -f (0oe)'^' h,,m', 



-V'p',i + IJLV'V^V^ = GijV'f TZij = v-p'{v'i,j + V'j,,) 

- 71, - e,j + V,, = G.fe^X + G,-fe^-f (<^z5,, + (t>nHij)e =v -^V^] 



[TZtj - Ttj - Stj)/v = gGik{v"v'k} + QGjk{v"v'k} + g{(pi5ij + <poHij)e Vi-j = -~v'{v[p',j +v'^p',i ) + iiv'{v'^V^v'- + v'^V^v'^) 
{TZ-T~e)/v''^gG,,{v':v';) + g{y, + cPoH)e 7^ = 7^.,/2; T = Tu/2; e = e,,/2; V = V,./2 

[Gi,'^]^ = 0; [G,fe^-<+G,fe<-<]^=0 <^.=Go5W; = ^[1 - ffW]; H = Hu/2 



hydrodynamics is so constructed that the effects 
of the mass flux and the density-spccific-volume 
covariance appear in closed form. Consequently, 
to represent these processes no explicit modeling 
is necessary, no additional equations need to be 
solved and the representation of these processes is 
mathematically and physically consistent. In VD 
moment closures these quantities require separate 
equations to be solved and their unclosed terms ap- 
proximated. 

3. Independence from the density model. The formula- 
tion for the momentum equation is independent of 
the particular functional form of the density PDF 
model, but requires the knowledge of the instan- 
taneous density, e.g. in the form of a SDE, such 
as discussed in Rcf. 2. This allows the momentum 
equation to be coupled to any mixing model, rep- 
resenting the fluid density. 

4. Consistent representation of the small scale ani- 
sotropy. Turbulence anisotropy is important in 
predicting the correct mixing state, mixing effi- 
ciency and the production of turbulent kinetic en- 
ergy. Compared to most second order moment clo- 
sures that model the large scale anisotropy and as- 
sume isotropic small scales, the developed stochas- 
tic equation represents anisotropy at both the large 
and the small scales. This is important in non- 
equilibrium, highly distorted, low- or modcratc- 
Reynolds-number or variable-density flows. 

5. Reduces to the original model for constant- density 



shear flows. Wc showed that the model coefficients 
can be designed so that the stochastic VD model 
reduces to the constant-density case in the limit of 
vanishing density fluctuations. 

Construction of a PDF model. A joint PDF model 
can now be constructed for variable-density turbulence 
where the stochastic density and velocity fields are cou- 
pled at the instantaneous level. The main ingredients 
for the hydrodynamics are collected in Table I, which is 
the essence of the paper. The VD momentum equation 
(first line) must be coupled to a fluid mass density PDF 
model, such as discussed in Rcf. 2, that represents con- 
servation of mass and the mixing. Then, in the velocity 
model, the tensors Gij and Hij , the mix metric function, 
g{9), and the kinetic energy dissipation rate, e, must be 
specified. Finally, the system must be augmented by an 
equation of state. The SDKs are discretized and numeri- 
cally integrated in time for a large number of Lagrangian 
particles (representing the fiow itself) in a Monte-Carlo 
fashion. Such a joint PDF model and its predictions for 
Raylcigh- Taylor flows are discussed in Ref. 21. 
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